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I. INTRODUCTION 


Random matrices play an important role in physics, engineering and mathematics men. While dealing with random 
matrices one often faces the problem of calculating eigenvalue spectra of sums and products of random matrices. This 
problem can be partially solved in the case of infinite invariant random matrices using methods of free probability 
which is a powerful technique developed recently hhiu. In practice, with some caution, one can use free probability 
results also as a good approximation for finite (large) matrices. 

Free probability calculus was introduced in operator algebra mu]. Later one has realized that it can be also 
applied to random matrices in calculations of moments of sums and products of large invariant hermitian random 
matrices m- In short, using free probability one can express the moment generating functions for a sum A + B or 
a product AB of invariant random matrices A and B in terms of the moment generating functions of these matrices. 
If a random matrix is hermitian the knowledge of its moments is often sufficient to reconstruct its eigenvalue density. 
Unfortunately the product P = AB , even for hermitian matrices, is generically non-hermitian. To get around this 
problem some authors |12] consider a related hermitian random matrix P' = \J~AB\fA, which has identical moments 
as P: Tr P' k = Tr P k . Unfortunately the eigenvalue density of P', except for a special class of matrices, is not the 
same as of P. So if one is interested in the eigenvalue density of P one has to either apply different methods H3IH3 
or to extend free probability calculus to the realm of non-hermitian matrices. For example, one might be interested 
in the eigenvalue density of the product of two independent GUE matrices. Clearly, the product is a non-hermitian 
matrix, its eigenvalues are complex and the eigenvalue density has a non-trivial support. Actually one can find the 
limiting eigenvalue density of this product to be p(z) = (27t| 2:|) _ ' 1 on the unit disk and zero outside the disk [T3l 1T4] . 

Non-hermitian random matrices have been studied in many contexts, including chaotic scattering |15lll8j . quantum 
chromodynamics with chemical potential |T9l{20], networks j23][24! , lagged correlations [2T[[22], matrix-valued random 
walk [25] ;25] and many others End]. As a result of these intensive studies one has developed analytic techniques to 
calculate the eigenvalue spectra of such matrices. In the heart of these techniques lies an electrostatic analogy mm- 
In this analogy eigenvalues of the random matrix play the role of positions of electric charges on two-dimensional 
plane and the Green’s function of electric field. Given the Green’s function one can calculate the eigenvalue density 
from the Gauss law [nidi did]. The Green’s function can be viewed as an subordinated element of a two-by-two 
matrix being an extended Green’s function. There were two independent approaches proposed for this extension: the 
first one was directly derived [29H521 from the regularization of the electrostatic potential mmm and the other 
one was obtained by a method of hermitization [55 -. I34j . Both the constructions are equivalent to each other and both 
can be formulated in terms of quaternions as one can see from the today’s perspective. 

The quaternionic nature of these equations was anticipated in [55] and worked out in [35, 30 where a very deep 
geometrical understanding of the regularization procedure was achieved. More precisely, in the calculations one 
usually introduces a regulator e. The modulus |e| can be interpreted as a distance from the complex plane. This 
means that calculations are in fact performed in the quaternionic space in the neighborhood of the complex plane and 
only towards the end of the calculations, when the limit e —* 0 is taken, the problem is projected back to the complex 
plane. The idea is similar to that used in the Sokhotski’s formula Q for one-dimensional delta function on the real 
axis except that now it is applied to two-dimensional delta function on the complex plane. In order to exploit this 
similarity it is convenient to recall the Cayley-Dickson construction. In this construction complex numbers are defined 
as pairs of real numbers and quaternions as pairs of complex numbers. In the one-dimensional Sokhotski’s formula 
one sends the second element of the Cayley-Dickson pair (imaginary part) to zero e —> 0 and regains one-dimensional 
delta function on the real axis while for the two-dimensional case one sends the second element of the Cayley-Dickson 
pair forming quaternion to zero and recovers the delta function on the complex plane [55, ,55 [ 15 0 . 

Yet before the importance of the underlying quaternionic structure of the extended resolvents was recognized 
[55] [55] one was able to derive the corresponding extended R transform (or equivalent blue function, sometimes used 
in physical literature) and to formulate the addition law for non-hermitian matrices [291433] . This law was then 
successfully applied to a number of problems, while the multiplication law for non-hermitian matrices was formulated 

much later 133 - 

Although quaternionic picture emerged already a decade ago [55] |55J [50 and has been applied to various problems 
[ 251 [55H5U] . it has not been commonly adopted. The reason for this is probably related to the fact, that the synthetic 
picture of what is the R transform for non-hermitian matrices was still missing. It was not clear how the R transform 
was related to the moments of the random matrix or what was the classification of the R transform in the simplest 
case of gaussian laws, etc. We fill up this gap here. 

We show that the quaternionic structure of the resolvent and the R transform is capable to incorporate all mixed 
moments of random matrix and its hermitian conjugate: ( j£TvX a X^ b X c ...), for an asymptotic class of random 
matrix models define by the probability measure of the type (44) for N — > oo. In contrast to hermitian case, though, 


these moments cannot be constructed from the moments of the eigenvalue density p( A) as they also depend on the left 
and right eigenvector correlations of the matrix X mm- So in a sense, the Green’s function and the R transform 
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encode also information about the eigenvector correlations. 

The link between free probability and random matrices was established through planar combinatorics m- In 
fact, the relation of random matrices to planar combinatorics was discovered long before free probability, by gaussian 
perturbation theory of invariant matrix models, which maps the problem of calculating moments of random matrices to 
enumeration of planar Feynman diagrams |43i I44j . Planar Feynman diagrams can be enumerated by field theoretical 
methods, in particular by Dyson-Schwinger equations, which relate various moment generating functions [441 [45]. 
These equations are very powerful, as one can use them to derive all relations between generating functions known 
from free probability. The idea is very universal and it can be also applied to non-hermitian matrices to derive 
relations between the quaternionic R transforms for independent random matrices A 1 B and the corresponding sums 
A + B [29H32] and products AB m- We follow this approach here. 

While presenting the material we first recall the hermitian case, which is well known and directly linked to free 
probability, and then we use it as a reference point while showing construction of the corresponding objects for 
non-hermitian matrices. We begin with the Green’s function and then turn to the R transform. Once we have 
recalled the definition of R transform for the non-hermitian matrices [35) 136], we will classify all possible forms of this 
transform for gaussian elliptic ensembles [45| . Then we extend the discussion to non-gaussian invariant measures (44) 
and link them to planar Feynman diagrams [Hi, 45]. We recall the addition law P2TJMTT21 and the multiplication law 
for non-hermitian matrices [37| and use them to evaluate eigenvalue density for a couple of examples of products of 
gaussian random matrices from elliptic ensembles. We shortly summarize the paper in the last section. In Appendix 
[A] we recall the Cayley-Dickson construction of quaternions. In Appendix [B] we discuss quaternionic representation 
of the two-dimensional delta function, in Appendix [C] we recall the corresponding representation of the sum of delta 
functions located at the eigenvalues of a given non-hermitian matrix m and finally in Appendix D we discuss planar 
Feynman diagrams [45H45] and recall how to use the diagrammatic technique to derive the basic relations ( 15|28 ) 
between the Green’s function and the R transform. 


II. QUATERNIONIC RESOLVENTS FOR NON-HERMITIAN MATRICES 


In this section we discuss Green’s functions (resolvents) for random matrix ensembles in the limit of large matrix 
size. We begin with the standard case of hermitian matrices and then generalize it to non-hermitian matrices. The 
generalization is closely related to the Cayley-Dickson construction of quaternions. We briefly recall this construction 
in Appendix [A] For hermitian matrices eigenvalues are real and the corresponding resolvents are complex, while for 
non-hermitian matrices eigenvalues are complex and the corresponding Green’s functions are quaternionic. 

Our goal is to calculate eigenvalue distribution of a random matrix H which is defined by a probability measure 
dp{H). Assume that H is hermitian H = H\ so it has real eigenvalues. In the limit of infinite matrix size N —> oo, 
the eigenvalue density of H is given by 


p(x) 


lim 

N—>oo 


l 

N 




2 = 1 



f 1 J v 

2 = 1 


(1) 


where A.;’s are eigenvalues of H. It is convenient to define the Green’s function (resolvent) which is a complex function 
of a complex argument z = x + ie 


G(z) = lim ( - H )- 1 \ = lim ( -J- . 

N— >oo \ N ! N— >oo \ N ^ Z — Xi 


N 1 


( 2 ) 


The symbol 1 stands for the identity matrix of size N. The eigenvalue density can be calculated from the resolvent 

© 


pi A) =-lim 3m G (A + ie) 

TV e->0+ 


(3) 


as follows from the Sokhotski’s representation of the one-dimensional Dirac delta function 

S(x) = -lim 3m-= — lim —- - . (4) 

7 r e->o+ x + ei tt e->o+ x z + e z 

The limit e —> 0 + projects the problem to real axis. All eigenvalues are located on this axis ©. 

One can now implement an analogous strategy for non-hermitian matrices. Eigenvalues are complex now, so one 
has to go beyond the complex plane. Using the Cayley-Dickson construction we first introduce a second auxiliary 
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complex variable, w, to define a quaternion q = (z, w) = z + wj (see Appendix|A| and a quaternionic Green’s function 
G(q). Once the Green’s function is determined one can send the second part of the quaternion to zero, w —> 0, and 
project it to the complex plane where eigenvalues reside. This is analogous to sending the imaginary part in the 
construction of eigenvalue density for hermitian matrices ©• „ 

Let us discuss this in detail. We use the notation of quaternions and of quaternionic matrices as in Appendix [A] 
We start with the quaternionic representation of the Dirac delta on the complex plane 


6™ (z) 


- lim SzF— 
7 r uj-vo oz z ■ 


1 


d 


- wj 


= — lim —F(z, w 

7T w-tO OZ 


I 1 = — lim 


M' 


7r w— >-0 


\z\ + \w 


\2 ■ 

2 ) 


(5) 


where F stands for the first Cayley-Dickson part of quaternion. It is analogous to the real part life of complex number. 
This representation is an counterpart of the one-dimensional delta function representation Q. The delta function is 
obtained by calculating the inverse quaternion close to the complex plane and eventually taking the limit by sending 
the second part of the quaternion to zero. The difference as compared to Eq. @ is that the imaginary part 3m 
is replaced by the first part F and that there is an additional derivative. In Appendix [B] we explicitly demonstrate 
that the expression <[5j) indeed yields a representation of the two-dimensional delta function. The representation ([5]) 
provides us with a guideline as to how to define the quaternionic Green’s function and how to calculate the eigenvalue 
distribution from it. 

The quaternionic resolvent for a non-hermitian matrix X is defined in a similar way as for hermitian matrices © 
mm- The Green’s function G(q) and its argument q are now quaternions [551 136]. Denote the first part of the 
Green’s function by G(q) and the second part by T(g), Q(q) = (G(q),T(q )) where q = (z,w ). In analogy to Eq. ([2]) 
we have 


Q (?) = (G(q), T(q)) = lim ( -^Tr b (zl - X, wt)~ 

N—too \ 


where Tr b is quaternionic trace defined in (A12). In the Pauli matrix representation (AlOl this equation reads 

S(q) = 


g(q) m \ _ lim n 

-T(q) G(q) J ~ N —+00 \ N b 


zl — X wl 
—wt zt — X t 


-i\ 


( 6 ) 


(7) 


In this representation the quaternionic trace Tr b is equivalent to taking trace of each of four N x N blocks separately. 
The matrix on the right hand side can be inverted 

zt — X wt \ 1 f (zt — X^)Hj^ —wFfft 1 


—wt zt — X t 


WFlr 


(zl - X)H^ 


where 


H l = (zt - X)(zt - A' f ) + \w\ 2 t 
H r = (zt - X^)(zt -X) + \w\ 2 t . 
We see that the first part of the quaternionic Green’s function is 


z 1-At 


G(q) = FG(q) = lim — (Tr(zt - X^HT 1 ) = lim — /Tr-- . lo 

n ^ ooN \ v ) l! n ^ 00 n\ ( z t-X)(zt-Xi) + \w\ 2 t 


and its derivative 
d 




= lim — < Tr 


N^oo N \ ((zt - X)(zt - Xt) + \w\ 2 t) ((zt - Xt)(zl - X) + |zc| 2 l) 


( 8 ) 


(9) 


( 10 ) 


( 11 ) 


In the limit of the second part of the quaternion q = ( z,w ) tending to zero, w —> 0, the expression in the brackets 
takes the form of a sum of Dirac deltas located at the eigenvalues of the matrix X 


N 


lim -^G((z,w)) = lim ( ^ ^ 5 (2) (z - A») } = p(z) 

0 7T OZ AT-S-oo \ Jy / 


( 12 ) 


so the eigenvalue density is regained in this limit similarly as in the one-dimensional case ([3]). The mechanism 
of localization of the expression in the brackets (11) at eigenvalues of the matrix X is analogous to that in the 
representation © of two-dimensional delta function. The details are given in Appendix [C] where we recall the 
argument given in m 
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III. QUATERNIONIC R TRANSFORM 

As in the previous section let us begin by recalling the hermitian case which is well known. It will serve as a 
reference point. The 1/z-expansion of the Green’s function § generates moments 


= lim ( —TV H r 
N—too \ N 


of the random matrix H : 


( 13 ) 


g(*) = £ 


m n 

7 n+l 


n —0 


The R transform is defined by the following equation [TO! 

G(z) = 


z-R(G(z)) ' 

The R transform is a generating function for planar connected moments K n of the matrix H 


( 14 ) 


( 15 ) 


^0) = E 


KnZ 


n —1 


( 16 ) 


The diagrammatic derivation [441 H51 of the relation between the Green’s function and the R transform is presented 
in Appendix D. The connected moments n n are also known as ’’free cumulants” or ’’non-crossing cumulants”. They 
are related to the moments m n as follows: = mi, k ,2 = m 2 — mj, K 3 = m 3 — 3 m 2 mi + 2 mf, K 4 = m .4 — 4 mini 3 — 

2 m\ + \Om\m 2 — 4 mf, etc. These relations can be derived by expanding (15). Note that the first three relations are 
identical as the corresponding relations between standard cumulants and moments known from classical probability. 
The reason why the R transform is important is that it is additive for sums of independent random matrices as 
discussed in the next section. 

Now consider a non-hermitian random matrix X. The situation is a bit more complicated since the 1/g-expansion 
of the quaternionic Green’s function ([ 6 ]) generates moments of differently ordered powers of the matrix X and its 
hermitian conjugate X'. 

It is convenient to rewrite the resolvent © as follows [JJ5J [5H] 




-1 


where 


A = 


A11 A12 
A 21 X 22 


X 0 

0 At 


and Q = q ® 1 where 1 is the N x N identity matrix and 


<1 = 


z w 
—w z 


( 17 ) 


( 18 ) 


( 19 ) 


is a matrix representing the quaternion q = (z, w ) (A7). We will index elements of the quaternionic matrix q by Greek 


letters q a p and consistently the positions of the blocks X a p in the matrix X (18) and in other block matrices. 
The 1 /q expansion of the Green’s function can be written as 


g(q) = lim / — Tr b ( Q - X)~ 

N—± 00 \ J\ 


= lim 

N—> 00 


N 


Tr b Q" 


+ lim 
N— KX> 


TibQ-^Q- 


lim 

N—>oo 


Tr h Q~ 1 XQ~ 1 XQ- 


( 20 ) 


or in the index notation as 


S{q)aC = q a l + E + E iToo ( 

M N ' P~/5e N ' 


+ ... . 


(21) 
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We can now define moments m^ of order n as multidimensional arrays with the following elements 


7771 n ^ 

'"'a 1 a2-..a2n-iC‘2n 


‘ ' ' '^» 2 n-l« 2 r 

JV -700 \ iv 


( 22 ) 


so we have 


G{q)a( = qj + TO ? 7 < 5 ^a /3 9 7 /+ • ■ ■ ■ 

/37 y07 5e 


(23) 


The last equation can be simplified since the off-diagonal blocks of the matrices X (18) are zero, so the only non-trivial 
moments are those involving only diagonal blocks: 


/^TrA’ aiai ... T Qria „ 


Combining that with Eq. (p3| we arrive at 


G{q) a <: — q a £ + qapQpc, + q aB q B ^ q~J +■■■ ■ 

0 0 j 


(24) 


(25) 


This equation is analogous to Eq. (14) except that now the moments are given by multidimensional arrays which 
encode information about all mixed moments. For example (M^) 12112 = limjv-»oo (^TrXXCY 2 X^). The first and 
the second moments are 


and 


Mj 1) = lim /— TrA), M 9 (1) = lim /—TrA' 1 ) 
1 TV—foo \N / 2 TV—7oo \N / 


= lim / AttAA ) , M$ = lim / AlVAA 1 > 


(26) 


H N—too \ N 


N—>oo \ N 


M 9 ( ? } = lim / —TrJYtX) , = lim / — TrA 1 A 1 ) 

21 jv-s-oo \ TV 7 ’ 22 " ■- x AT 7 




N—t oo \ AT 


(27) 




respectively. Generally the array representing the n-th moment has 2” elements, but not all of them are independent. 
For example M-fy 1 = M 22 \ M 2 ^ 


The quaternionic R transform is defined in the same way as the R transform for hermitian matrices (15) [2911301 


G(q) = 


q-n(G (q)) 


(28) 


except that it is a quaternionic function. It generates planar cumulants K which are now multidimensional arrays 
in contrast to the hermitian case where they are real numbers (16) 


^ (uOaC — + ^aC K a 0£qa0q0(, + K a 0^ ( .q ol 0q0r y q 1 ^ .... 

0 07 


The first cumulant is 


and the second one is 


k[ 1] = lim /— TrA'V K { 2 ] = lim /—TrA'A 
1 TV-foo \ N / 2 N^oo \N / 


(29) 


(30) 


< = - xi ^ x - *»>) • - V 

Ag> = »lTo„(v li(x - : ' 1),(x - 3:1) )' = -*»>’) . 

We used here a shorthand notation x = K\ 11 for brevity. 


(31) 








7 


IV. GAUSSIAN ELLIPTIC ENSEMBLES 


In this section we consider the class of gaussian elliptic random matrices whic h a re defined by the condition that 
the third and higher cumulants vanish: K^ = A'U) = ... = 0. The R transform (29) is a linear function in this case 


Mq)ac = Ki a )s °c + K$qc'<; ■ 

The last equation can be rewritten as a matrix equation 

K(q) = I< {1) + A (2) o q , 


(32) 


(33) 


where the symbol o denotes the Hadamard product and , K^ are two-by-two matrices representing the first two 
cumulants 


K& = 




I< (2) = a 2 


He 2irt> 1 \ 

1 ) 


(34) 


The first cumulant is parametrized by two real parameters: 3 ftex and Qmx while the second one by three: cr, fi and 
</>. The eigenvalue density of the corresponding random matrix is uniform on an ellipse on the complex plane. The 
parameter x is the position of the ellipse center, a 2 = af + a\ is the sum of squares of lengths of semi-axes, fi is an 
eccentricity parameter defined as /i = (a 2 — erf)/(cr 2 + erf), and is the angle between the first semi-axis and the real 
positive semi-axis. We define a family of standardized elliptic laws by setting x = 0 and a = 1 and <j> = 0. All elliptic 
laws can be obtained from the standardized one by shifting it by x, rescaling by a factor a and rotating by <j>. The 
standardized family is parametrized by one parameter /i. The R transform is 


n 


z w 
—w z 


fl 1 

1 n 


z w 
—w z 


[XZ w 
—w nz 


Using the notation with quaternionic units (A3) we can write this equation as 


(35) 


lZ(z + wj) = + wj . 


(36) 


We see that the R transform is given by a very simple linear map which multiplies the first Cayley-Dickson part (A6) 
of the quaternion by H'- A(7?.(g)) = /i-F(g) and leaves the second part intact S(TZ(q)) = S(q). For hermitian matrices 
= X, h = 1 and 7 Z(q) = q. For anti-hermitian matrices X t = —A', h = — 1 and 7 Z(z + wj) = —z + wj. For Ginibre 
matrices m H = 0 and lZ(z + wj) = wj. The R transform for a general elliptic law reads 


lZ(z + wj) = x + a 2 (/re 2 *^ + wj) 


(37) 


One should note that the rotation acts only on the first part of the quaternion, while the rescaling on both. 

One can easily derive the probability measure for elliptic gaussian random matrices. Elliptic random matrices can 
be constructed from gaussian hermitian matrices. Let Hi and H 2 be two independent N x N gaussian hermitian 
matrices defined by the following product probability measure 

dH{H u H 2 ) (x D7L 1 D7J 2 e-^^'K + ^ 2 ) , (38) 

where DH\ and DH 2 are flat measures for hermitian matrices each having N 2 real degrees of freedom. We skipped 
the normalization constant in the last equation. As the measure factorizes, we have 


^^Trf7 a ^ = 0 , ^ ^TrH a H b ^ = 5 ab . (39) 

for a = 1,2, b = 1,2. Let us define A as a linear combination of these matrices 

X = xt + e*+ (tn-ffi + ia 2 H 2 ) . (40) 

This combination is generically non-hermitian and has 2 N 2 degrees of freedom. Its hermitian conjugate is 


X^ = xl + e (criHi — itj 2 H 2 ) ■ 


( 41 ) 





It is easy to calculate the one-point and two-point correlation functions (30) and (31) for A and X'. One gets the R 


transform (34) or equivalently (37). In order to obtain the measure for X one can invert the two last equations for 
Hi and H 2 and insert the result to (|38|). Without loss of generality we set x = 0. We have 


Hi = 


e-^X ■ 


- e+^A't 


2(71 


ill = 


e~ ill, X - e+^Xi 


2 ia 2 


and 


dn 0 (X) ex DX exp ^ ^ 2) Tr (XX* - | (e" 2 ^A 2 + e + 2 i *A t2 ))^ 


(42) 


(43) 


where DX is a flat measure for non-hermitian matrices which has 2N 2 real degrees of freedom. For x 0 the matrix 
X should be replaced by X — xl. It is the most general form of the gaussian elliptic measure. For any /r £ (—1,1) 
the distribution is uniform on an ellipse. The limit fj, —> ±1 should be taken carefully. In this limit the exponent in 
(43) changes to a delta function 5(e _l ^A e +l ^X') which is responsible for a reduction of the number of degrees 

of freedom from 2 N 2 to N 2 . The width of the ellipse shrinks to zero and the resulting eigenvalue density has a one 
dimensional support being an interval on the complex plane. The eigenvalue density is not uniform on this interval 
but it is given by the Wigner semicircle law. It is worth noting, that those are only two limiting cases, while density 
experiences an interesting crossover regime when fi ~ ±1 -k |4&) . 


V. NON-GAUSSIAN ENSEMBLES 

In this section we go beyond the gaussian regime. Consider random matrices defined by the following invariant 
probability measure 

dfi(X) cx (44) 

where V(A, A't) is a two-matrix polynomial such that TrV(X, At) is real and bounded from below. This means 
that there exist a real number r such that TrR(A, A') > r for any matrix A'. Such measure is invariant w.r.t. 
transformation A —> UXU ', but such invariance is not sufficient to express the measure only in terms of eigenvalues 
like in lrermitian case. The left and right eigenvectors are not simply complex conjugates of each other, but their 
correlations carry an non-trivial information. The polynomial V(X, A') is constructed as a sum of terms of powers 
of A and A^ 


gX kl X* jl ...X k "XV” 


(45) 


with complex coefficients g called coupling constants. Different terms may have different coupling constants. Terms 
which can be obtained from each other by a cyclic permutation of matrices have identical trace and are thus equivalent 
from the point of view of the measure (44). It is therefore convenient to divide all terms into equivalence classes of 
products of A and A^ which can be obtained from each other by a cyclic permutation. For example the following 
products A 2 A 3 ’!'A 2 , A 4 A 3 ’f, X^X 4 X 2 ^ belong to the same equivalence class and have the same trace TrA 2 A 3 'l'A 2 = 
TrA 4 A 3 l = TrAfA 4 A 2 l. We call these equivalence classes periodic chains. We also define dual chains. A chain 
is said to be dual to a given chain if it is obtained by swapping positions of A’s and A^’s in the chain. The 
requirement that the trace of the potential is a real number means that the potential V(A, A^) must contain pairs 
of dual terms with complex conjugate coupling constants. For example the terms gX 2 X 1 and gX^ 2 X are mutually 
dual with proper coupling constants. For self-dual chains, as for instance <?A 2 A 2 ^, we have g = g, therefore those 
coupling constants must be real. For example the most general form of the third order terms in the potential is 
giX 3 + <71 A'i ' 3 +g 2 X 2 X^ + g 2 A^ 2 A with two complex coupling constants denoted here by gi and g 2 . The fourth order 
potential has two dual pairs 33A 4 + <73A'i' 4 and g±X 3 Xi + g 4 Xi 3 X : as well as two distinct self-dual terms g$X 2 X t 2 
and g§XX'XX' with real coupling constants g$ and g§. 

For convenience we assume, without loss of generality, that the measure is centered, that is (ATrA) = 0, so we 
can skip linear terms in the potential. Next step is to split the measure into the gaussian part (43) and to treat the 
residual part as a perturbation of the gaussian measure [44] 


dfi(X) = d» 0 (X)e- NTrV «( x ’ x ^ = dy. 0 (A) - TVTrV^A, At) + iA 2 Tr 2 VR(A, A*) + ..}j 


(46) 


where dg, o(A) is the elliptic gaussian measure (43) and Vr is the remaining part of the potential which contains the 


third or higher order terms. In this way, the average of an observable 0{ A) is calculated as an average over the 
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FIG. 1. Three possible propagators coming from gaussian part dfi o (A') of the probability measure. 


gaussian measure d/io(A) but of a new observable modified by terms coming from the perturbative expansion of the 
residual part 


(O(X)) = (o(X) (l - NTrVn(X, AT+) + * N 2 Tr 2 V R (X, X*) + ■••)) 


(47) 


Now the calculation is reduced to gaussian integration. Thanks to the Wick’s theorem it can be mapped onto the 
problem of Feynman diagram enumeration which simplifies in the limit N —> oo to enumeration of planar diagrams 
[43l [44]. We will not give the details here, the interested reader can find them in m ■ We only give the Feynman 
rules for this model. Diagrams are constructed by drawing lines between vertices. The lines, called propagators, are 
deduced from the gaussian part d/io{X) ( [43] ) . There are three different cases shown in Fig. [T] which correspond to 
two-point functions (XX)o, (XI*)o (XbXQo. The vertices come from the expansion of the residual part. They 
correspond to different periodic chains. As an example, we show in Fig. [2] a diagram consisting two third order 
vertices contributing to (XXX' X'~). For each such diagram, there will be a dual one obtained by replacing matrices, 


X o X\ and vertices, f> jj, to their dual ones. The sum of all contributions from connected diagrams with n 
external lines is equal to planar cumulant of order n. Sometimes it is called a non-crossing cumulant or free cumulant 
of order n. One can symbolically denote cumulants as connected averages ((-^TrX a X 1 ' b A' c ...)). A cumulant of order 
n is an array which contains connected averages for all distinct periodic chains of length n = a + b + c+ .... The R 
transform is a generating function of such cumulants (291. The main relation (151 between the Green function and R 
transform can be derived from the Dyson-Schwinger equations for planar diagrams. The details can be found in m- 
The R transform has an important property. The R transform of a sum A + B of two independent random matrices, 
each defined by a probability measure of the type (441, is additive 


Xa+b{<i) = X,A(q) +77s(g) . (48) 

This property can be easily understood in terms of Feynman diagrams. The measure d/j,(A, B) = d[i{A)dfi{B) 
factorizes for independent matrices and also its gaussian part does d^o{A, B) = dfj,o(A)dfio(B). As a consequence, 
the gaussian part does not contain mixed AR-terms and therefore the mixed propagators vanish ( AB) 0 = (AB^) 0 = 
... = 0. This means that there are no lines in the Feynman diagrams which would directly connect vertices of type 
A and B. In effect, any connected diagram contains either vertices solely of type A or vertices solely of type B. In 
the former case they come from TZa(q) and in the latter one from 7 Zs{q)- Thus the generating function for connected 
diagrams TZa+b{q) splits into separate contributions from diagrams of type A and of type B: TZa (<?) + TZb{q)- This 
is the standard diagrammatic interpretation of the additivity of cumulant generating functions. 

One can work out consequences of the absence of mixed propagators also for the product AB of independent 
matrices. We do not show this calculation here and refer the interested reader to m • The multiplication law can be 
written in the form of three equations. 


TZab {Gab) = [Xa {Gb)] L [TZb {Ga)] R , 

[Ga] R = Gab [TZa ( Gb)] L , (49) 

[Gb] L = [TZb {Ga)] R Gab, 


Since now the R transforms are quaternionic (or matrix-valued) the order of multiplication matters. The quaternions 
Ga = Ga{z ), Gb = Gb(z), Gab = Gab(z) in the last equation correspond to the Green’s functions for A, B and 
AB projected to the complex plane, that is calculated for quaternion z = (z, 0) having the second part equal zero. 
Quaternions tagged in the last equation by L or R correspond to left or right rotated copies of quaternions, given by 
the following transformations: 

[Q] l = UQU\ 

[Q] r = U^QU. 

where U = diag (e*^/ 4 , e ~ I?i / 4 ) is a unitary matrix constructed from the phase of the complex number z = re 1 ^. 


(50) 
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FIG. 2. Example of a diagram consisting of all four kinds of third order vertices (light gray circles), contributing to cumulant 
of order 4: ((TrX 2 X' 2 )). The big dashed circle emphasises symmetry with respect to cyclic permutations under trace. 


For hermitian matrices the quaternionic set of Eqs. 
and Rab being complex functions 


l uy i smipnnes 






Rab {z) = Ra (w) R b (v) , 

v = zR a (w) , (51) 

w = zRb (v) . 

Moreover, if additionally the first cumulants of A and B are non-zero the last set of equations can be concisely written 
in this case in terms of the S transform [0] 


Sab(z) = S A (z)S B (z) (52) 

which is related to the R transform as follows m 

Sa{z) = Ra{zS a {z)) ° r Ra{z) = Sa(zR a (z )) ' (53) 

Unfortunately, there is no corresponding S transform representation of the multiplication law for non-hermitian 
matrices. 

We finish this section by giving the transformation law for the R transform under multiplication of the matrix by a 
complex number: A —> A' = aA. Let us denote by d a quaternion whose first Cayley-Dickson part is equal a and the 
second part is zero: a = (a,0). One can see [351136] directly from the definition of the quaternionic Green’s function 
([6]) that 


g aA (q) = Ga(& 1 q)a 1 . 


( 54 ) 
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FIG. 3. Schematic representation of addition law in case of sum of Ginibre and GUE matrices, which results in Elliptic matrix. 


If follows from Eq. <[15j) that the R transform obeys the following transformation law under multiplication of the 

n aA {q) = all A (qa) ■ 


matrix by a complex number 


Additionally if we include constant shifts in the transformation law A —► A' = xt + aA we obtain 

H A '{q) = x + oill A (qa) . 


(55) 


(56) 


In particular when we apply this transformation to the standardized elliptic gaussian law (36 1 we obtain the general 
form (l37|). 


As an example of the application of this transformation law let us derive the R transform for X = -U (A + iB) 
where A and B are standardized independent gaussian hermitian matrices: lZ A (q) = IZs(q) = q- Using the addition 
law (48) and the transformation (l55| for rescaling we easily find in agreement with [23 [3B] 


K x ((z,w)) = i( z,w ) + i(i, 0 )( 2 ,w)(i, 0 ) = ( 0 ,w) , 


(57) 


where the argument of the R transform is written as a Cayley-Dickson pair q = (z,w). In the last step of the 
calculation we applied Eq. ( |A5[ ). As expected, we reconstructed the R transform for the Ginibre matrices. For 
pedagogical reasons we rewrite the last equation in the matrix representation (A7) which is better known 


Tlx 


z w 
—w z 


1 f z w 

2 V —w z 


+ l( i 0 

+ 2 l 0 -i 


2 w 
—w z 


i 0 
0 -i 


0 w 

—w 0 


(58) 


VI. APPLICATIONS OF QUATERNIONIC R TRANSFORM 


In this section, as an illustration, we give several examples of how to calculate eigenvalue densities of sums and 
products of gaussian elliptic random matrices using the quaternionic R transform. 


A. Sums of gaussian random matrices 


It is easy to see that a sum C = A + B of gaussian elliptic matrices A and B is again an elliptic gaussian random 
matrix. Denote the R transforms (371 for A and B as 1Z A {z + wj) = x A + a\ ( y ^ A e 2l ^ A z + wj ) and 7 Ib(z + wj) = 
z + wj ). As follows from the addition law (48) the R transform for C has exactly the same form 
(/Jce 2 ^ 0 z + wj) with xc = x a +xb, = a A + a 


xb + (uBe 21 ^ 13 


Tl c (z 


vj) =xc + <?c 


I, nee** = {n A (J 2 A e i ‘t >A +HBa 2 B e i * B )/{a 2 A + 


<Jb)- In other words, the elliptic gaussian laws are stable under addition. In particular if we add a Ginibre matrix X 
(x A = 0 ,fi A = 0 ,<j a = 1 ,<j> A = 0) to a hermitian GUE matrix H [1( (x A = 0 ,fi B = 1, cfb = 1 ,4>b = 0) we obtain a 
matrix C = X + H with the R transform 7 lc(z + wj) = z + 2 wj. The corresponding eigenvalue density is uniform 
on an ellipse located at the origin of the complex plane with the longer semi-axis of length %/2 on the real axis and 
the shorter one of length 1 on the imaginary axis. This is schematically shown in Fig. [3] 


B. Products of gaussian random matrices 


The eigenvalue density of a product of gaussian elliptic random matrices can be calculated from the general multi¬ 
plication law Eqs. (49) and the Eq. (28 1 


Gab {z) 


1 

z - 7 l AB ( G A b ( 5 )) 


(59) 
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which relates the Green’s function and the R transform for the product AB. From these equations we obtain the 
Green’s function Gab ( z) on the complex plane that is for quaternions having the second part equal zero z = (z,0). 
Then we extract the first part of the quaternionic Green’s function 

G ab (z) = FG ab (z) (60) 

which is a complex function G A b(z) and eventually we use Eq. ( fl2| to calculate the eigenvalue density 

ldG AB (z) 

PAB{z) = ;—a^- (61) 

Let us give a couple of examples. For convenience we introduce a shorthand notation. We denote the standardized 
elliptic matrix with the eccentricity parameter /i by E(/i) 1461 . The matrix E(/i) has the following quaternionic R 
transform TZe(h){z + wj) = /j, + wj , or in the matrix representation 


n E(v) 


z w 
—w z 


{IZ w 
—w flZ 


(62) 


We reserve a special notation for the standardized Ginibre matrices X = E( 0) and the standardized hermitian GUE 
matrix by H = E( 1). Let us note that the Ginibre ensemble is invariant under rotation X —► X' = e^X. Indeed, 
using the transformation law (55) we see that 




z w 

—w z 


0 \ ( 0 w \ ( e** 

0 e-** ) \~w 0 J \ 0 



(63) 


the R transform stays intact unter the multiplication of the Ginibre matrix by the phase factor e 1 ^. Now we are ready 
to consider examples. 

We use the Pauli matrix representation for practical calculations, because complex algebra is familiar to all readers 
and is well implemented in many software packages, which is not the case for quaternions. 

The first example is the product of A = al + bX i and B = ct + dX 2 

AB = (al + bXx) (cl + dX 2 ), (64) 


where a,b,c,d € C, b, d ^ 0 and X\, X 2 are independent standardized Ginibre matrices. The question we ask here is: 
what is the eigenvalue spectrum of the product? For such a pr oduct one can easily see, using the invariance of the 
Ginibre ensemble under the multiplication X —► X' = e l ^X (63), that the following redefinition of parameters: 


s = 


t = 


u = |M|e iArg(ac) , 


(65) 


reduces the calculations to the problem of finding the eigenvalue density of a matrix 


- (si + (tl + X2 ), 


( 66 ) 


where s,t £ R+. The parameter u scales the eigenvalue density with \u\ and rotates it around the center of the 
complex plane by Arg (u) angle. This means that it is sufficient to consider products of the type 


AB = (si + X ± ) {tt + X 2 ) 

without loss of generality. Using the addition law ([48]) for matrices A = si + X and B = tt + X we have 


and 


TZa 


Kb 


W B 

VB 

-V B 

w b 

W A 

VA 

-VA 

W A 


s v B 
~VB S 


t VA 
-V A t 


(67) 


( 68 ) 


(69) 


respectively. Sinc e s and t are real we may omit their complex conjugation. Inserting Eqs. ( 68|69 ) to the multiplication 
law formulas (49) and writing 2 in the polar form 2 = re 1 ^ we get the R transform for the product: 


Kab (Gab) = 


st — v a vb£ 


i(f> 


tv A e 


—i(p/2 


SVbG 


40/2 


-tv A e irl,/2 - sv B e~ i ^ / ' 1 st - v A v B e"^ 


(70) 
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Now using Eq. (591 we obtain the quaternionic Green’s function on the complex plane. Its first Cayley-Dickson part 
reads 


Gab(z) = - 


st — re l<f> — vavb 


s 2 t 2 


■t 2 \v A y 


\V B \ 


'■{vavb 


■ vavb) + \va\ 2 \vb\ 2 - 2str cos 4> 


(71) 


Now we have to find va and v B as a function of complex argument z (or equivalently r and <j>). To this end we utilize 
the second and third quaternionic equations in (49) that are equivalent to the following four independent complex 
equations 


t\v A \ + e®’’ (w A r + v A v B (s + w A )) = t( 1 + sw A ), 
s\v B \ 2 + e l4> ( w B r + VAVB (t + w B )) = s(l + tw B ), 
2 ^ sw A v B + V B \v A \ 2 + v A r - Vb') = tv A (s + W A ) , 
~ l(f> (^-tw B VA + VA |ub | 2 + v B r - VaJ = SVB (t + W B ) ■ 


(72) 

(73) 

(74) 

(75) 


The last two equations have a trivial solution: v A = v B = 0. It corresponds to the holomorphic part of Green’s 
function, G (z) = (z — st ) _1 , valid outside the eigenvalue domain. 

By some algebraic manipulations in Eq. (741 we may express v B as 


VB=VA IAAF 


(76) 


which means that combinations v A v B and v A v B do not depend on phases of va and v B - Indeed an inspection of 
the expression for the Green’s function ( 7T|) sh ows, that it depends on va and v B only through the moduli IiuU'WbI- 
Their phases are also irrelevant in Eqs. ( 74|75 ). In fact, this is true for any product of elliptic ensembles. For the rest 
of the discussion those phases will be set to 0, so v A = v A , v B = v B and v A ,v B £ IR + . Solving Eqs. ( |72|73|74|75[ ) for 
va and v B yields: 


v A (— 1 +v\ + s 2 )(v b + t 2 ) + (-1 + 2 v\)v B r + v A r 2 - 2 stv A r cos <j> = 0 , 
i > b (—1 + v% + t 2 ){v\ + s 2 ) + (-1 + 2 v%)v A r + v B r 2 - 2stv B rco$,(j) = 0 . 


(77) 

(78) 


Note that those equations are invariant under the simultaneous change of sings va —t —va and v B —t — v B which is 
a remnant of the symmetry of v’s mentioned above. In effect if va and v B is a solution to these equations, then also 
v’ A = — va and v’ B = — v B is a solution that gives the same Green’s function, leading to the same eigenvalue density. 

The eigenvalue density is represented as a function with a compact support on the complex plane. Let us focus 
on the outline (contour) of this support. We may calculate it by matching the holomorphic (va = v B = 0) and 
non-holomorphic solutions. To do so, we solve quadratic Eq. for v B , and pick the leading term of the solution, 
that vanishes in the limit v A 0. Then by inserting this solution into Eq. and performing the limit va —t 0 we 
are left with the following equation 


( s 2 t 2 — s 4 t 2 — s 2 t 4 + s 4 f 4 ) + (2 s 3 t cos (f) + 2 st 3 cos </> — 4s 3 f 3 cos (j))r+ 
+( — 1 — s 2 — t 2 + 2 s 2 t 2 + 4 s 2 t 2 cos 2 (j>)r 2 — 4 st cos (j) r 3 + r 4 = 0 , 


(79) 


which may be treated as a fourth order polynomial equation for r and solved for any value of <f> giving the contour of 
the eigenvalue distribution. 

Now the remaining part is to numerically solve real polynomial Eqs. ( f77|78l ), e.g. on a lattice inside the previously 
calculated domain, plug the results into Eq. ( pflj ) and calculate eigenvalue density according to Eq. (12). The results 
for several different values of s and t are shown in Fig. [4j together with comparison to numerical simulations. The 
agreement is very good. Due to the finite size effects the edges of the support are smeared and instead of a sharp 
edge on the contour line there is a crossover region of width that decreases with N where the eigenvalue density falls 
off to zero. 


The equation (79) is general in the sense that it holds for any t and s. The special symmetric case t = s can be 
solved in a simpler way. The symmetry A <—> B which holds in this case significantly reduces the complexity of 
the problem because one may at the beginning of the calculation set w A = w B and va = v B and reduce the number 
of unknown variables. Two special cases: t = s = 0 and t = s = 1 of this symmetric problem have been discussed 
previously m and can be used as a cross-check of our general formula ( |79| . The case t = s = 0 corresponds to the 
product of two standardized Ginibre matrices with the contour given by the equation m 


r z - 1 = 0 . 


( 80 ) 
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FIG. 4. (Color online) Numerical eigenvalue distributions for matrices (sl + X)(tl + X) 100 x 100 for (a) s = 0.50 ,t = 0.75 (c) 
s = 0.9,t = 1.2 and (e) s = 1.2,t = 1.3 compared with theoretically evaluated ones in (b),(d) and (f) respectively. The edge of 
the distribution is shown on the plots (a),(c) and (e) (red curve). Each histogram is made from 10' eigenvalues. 


The case t = s = 1 corresponds to the product of two standardized Ginibre matrices shfited by unit matrix and has 
the contour given by the equation 

r 2 — 4 cos (j> r + 4 cos 2 <j) — 1 = (r — 2 cos <j> — l)(r — 2 cos <j> + 1) = 0. (81) 

The first one gives a circle and the second one a Pascal limagon, in agreement with mm- 
The second example we consider is a product 

AB = (1 + £d(/r))(l + E 2 (h)), (82) 

of two independent identically distributed shifted elliptic matrices A = 1 + and B = 1 + -E^q,). We follow 
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exactly the same procedure as in the first example except that now we exploit the invariance of the problem under 
the exchange Ei(/i) <—> E-z(n) to make a symmetric Ansatz wa = w B = w and v a = vb = v at the beginning of the 
calculations. This simplifies the problem remarkably. The R transform for a shifted elliptic matrix is given by (371 


n 


W V 

—v w 


1 + W/i V 
— V 1 + W/i 


and for the product: 

TIab (Gab) = 

The complex Green’s function is therefore equal 
Gab (z) = 


(1 + w/i) 2 - M 

-v (e-^/ 2 (l + w/i) + e^/ 2 (l + w/i)) 


v (e "^/ 2 (1 + wn) + e ^/ 2 (1 + w/i)) 
(1 + w/i) 2 - \v\ 2 e-^ 


(83) 


(84) 


z i(t> (1 + w/i) 2 -r- 


(:r — e 1 ^ — w 2 e l ^/i 2 — 2 e^wr) ^(1 + w/i) 2 — re ®^] — |u| 2 e 1 ^ 
The four complex equations for w and v reduce to two 


^ — |u | 2 e®^ ^2 ^1 + r + fi , ^ 


1 + r + fi [w + w + fi \ w\ 


■)) + w) 


o ir t> 


(wr + (1 + w + w/i) |«| 2 ) = (1 + w/i) ^1 + w + w 2 /i — |u| 2 ^ , 
g z 0 + r + |u | 2 — w (1 + w/j)^J = v (1 + w/i) (1 + w/i + w). 


(85) 

( 86 ) 

(87) 


v = 0 is again a trivial solution. A non-trivial solution is obtained by the procedure of matching holomorphic to 
non-holomorphic solutions by linear approximation in v in the very same way as in the previous example. We find 
equations describing the boundary of the eigenvalue distribution 


re l ^w = (1 + w/i) (l + w + w 2 /i) , 
e®^ (—1 + r — w (1 + w/i)) = (1 + w/i) (1 + w/i + w), 


( 88 ) 

(89) 


which may be solved for any value of <j> € [0,27t) giving the desired contour. Then one may numerically solve the 
polynomial Eqs. ( 86|87 1 inside the eigenvalue domain. The results for different values of /i, are presented in Fig. [5] 
supported by numerical simulation of random matrices. Theoretical predictions match the simulation up to finite size 
effects. 

Finally the last example to be analyzed here is the product of the form: 

AB = (1 + H)(l + A), (90) 

of a shifted standardized hermitian GUE matrix A = t + H and a shifted standardized Ginibre matrix B = t + X. 


The R transforms are special cases of the elliptic one (831, for /i = 1 and /i = 0 respectively. The calculation procedure 

(91) 


TIab(Gab) = \ _ iH2( 


is analogous to the previous ones. The R transform for the product is: 

1 + wb - VAVBe l<t> e - ®^/ 2 (va + wbva + nge®^) 

-e^' 2 (va + wbva + uge - ®^) 1 + w B - VAVse~ l<l> 

The Green’s function reads 
Gab (z) = 

1 - re~ i4, +w B ~ e-^VAVB 

2Re (re®^ (1 + w B )) - ( r 2 + ^2Re (w B ) + |wb| 2 ) (l + lu^l 2 ) + 2rRe (vavb) + (l + |^| 2 ) (l + \vb\ 


(92) 

(93) 


The solution va = v B = 0 again gives a holomorphic Green’s function. Then, in similar way to the previous two cases 
it is sufficient to consider linear terms to obtain eigenvalue domain. One is left with three independent equations 
relating wa, wb, r and <j> 


(1 + w b +w 2 b ) e-** 


w B 

wa = w B (1 + w B ), 


0 = -1 + \w B \ z + w B e 2 ®^ (1+wb+ w B ) ■ 


(94) 

(95) 

(96) 
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FIG. 5. (Color online) Numerical eigenvalue distribution for matrices (1 + i?i)(l + E 2 ) 100 x 100 for (a )fj, = 1/3 (c) ^ = 4/5 
and (e) fx = 23/25 compared with theoretically evaluated ones in (b), (d) and (f) respectively. The edge of the distribution is 
shown on the plots (a), (c) and (e) (red curve). Each histogram is made from 10' eigenvalues. 


It is pointless to give explicit solutions for w’s, as they are given by long and cumbersome formulas. One can 
nevertheless implement them to a numerical code to evaluate the position r of the density border for any (f) £ [0, 27t) 
through Eq. (94). It turns out that the solution has two branches, each forming a closed loop around a different 
portion of the eigenvalue density. The loops touch at a single point at the origin z = 0. The result and the comparison 
to numerical simulations for finite matrices is presented in Fig. [6] The finite N results based on the MC simulation 
converge to the theoretically evaluated domain while the size of matrix is being increased. 
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FIG. 6. (Color online) Numerical eigenvalue distribution for matrices (1 + H)( 1 + X ) of size (a) 50 x 50 (b) 100 x 100 and (c) 
200 x 200 compared with theoretical prediction for the edge of the distribution (red points). Each histogram is made from 10' 
eigenvalues. The eigenvalue domain is determined by closed regions given by analytic solutions of Eqs. (194} and (961. 


VII. CONCLUSIONS 


One can systematically define the Green’s function and the R transform for non-hermitian random matrices given 
by invariant measures of the form (44) and calculate the corresponding eigenvalue densities in the limit N —>■ oo. The 
Green’s function and the R transform have an analogous structure to those for hermitian matrices known from the 
standard free probability calculus, but they are quaternionic. We demonstrated how to derive these functions using 
the Cayley-Dickson construction. 


The addition and multiplication laws for invariant non-hermitian matrices can be deduced from Feynman pertur¬ 
bation theory for matrix models with the invariant measures. These laws provide us with a practical method to 
calculate the limiting eigenvalue densities for sums and products of non-hermitian invariant random matrices in the 
limit N —^ oo. 


If one combines the method discussed in this paper with the trick of linearization |25l 149] one can also apply it to 
determine eigenvalue densities of more complicated matrix polynomials. So far this method could only be applied to 
polynomials which were hermitian by construction as for instance AB + BA + A + B or ABC + BCA + CAB for A, B 
and C being hermitian. This required a special adaptation of the linearization method ESI nu. Now, this constraint 
can be removed and one can treat non-hermitian polynomials as well, as for instance AB + A + B or ABC + BCA. 

As mentioned before, the non-hermitian products of random matrices appear in multiple contexts including iterative 
stochastic evolution of linear systems [52H54] , capacity of complicated MIMO networks [55] [55] or analysis of multi¬ 
dimensional data with asymmetric correlations, like e.g. time series of stock prices ED ED HD HHi and others. Many 
of such problems were tackled with methods of free probability in the hermitian case. Now, by the use of the 
quaternionic R-transform, they can be extended towards the non-hermitian regime which is often computationally 
harder, but physically better motivated. 
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Appendix A 

In this Appendix we recall the Cayley-Dickson construction of quaternions and discuss an extension of this con¬ 
struction to matrices. 

Quaternions q are linear combinations 


q = a + bi + cj + dk (Al) 

of quaternionic units 1 ,i,j,k with real coefficients a, b, c, d. The quaternionic units form the Hamilton basis which is 
defined by the the following multiplication rules 


i 2 = j 2 = k 2 = ijk = —1 . (A2) 

Quaternions can be represented as ordered pairs (z, w) of complex numbers z = a + bi, w = c + di 

q = a + bi + (c + di)j = z + wj = (.z, w) . (A3) 

The conjugate of q is q* = a — bi — cj — dk = z — wj = (z, —w). Quaternion addition is defined by 

(z, w) + (v, y) = {z + v,w + y) , (A4) 


and multiplication by 


(z, w)(v, y) = (zv — wy, zy + wv) . (A5) 

It is useful to define auxiliary functions F and S which allow one to extract the first and the second element of the 
Cayley-Dickson pair q = ( z, w) 


z = F(q) , w = S(q) . (A 6 ) 

These functions play an analogous role for quaternions as the real and imaginary parts for complex numbers. We 
refer to them as to ’’first part” and ’’second part” of quaternion. Quaternions have a matrix representation in terms 
of Pauli matrices <r,. A quaternion q = a + bi + cj + dk= (z, w) can be mapped to two-by-two complex matrix (we 
use the same letter to denote this matrix) by 


q = aa 0 + b(icr 3 ) + c(*<r 2 ) + d(ia\) 


a + ib c + id 
—c + id a — ib 



(A7) 


The first row of the quaternion matrix representation (A7) can be identified with the Cayley-Dickson pair ( z,w ). It 
is easy to check that the matrices ao,ia 3 ,ia 2 ,i<Ji form the Hamilton basis (A2) 


(zct 3 ) 2 = (*cr 2 ) 2 = ( ivi ) 2 = (io- 3 ){ia 2 )(icri) = -cr 0 . 


(AS) 


Quaternion addition corresponds in this representation to matrix addition while quaternion multiplication to matrix 
multiplication: 


2 w 
-w z 


v y 
-y v 


zv — wy zy + wv 
—zy — wv zv — uiy 


(A9) 


as c an be seen by comparing the upper rows of matrices in the last equation with the Cayley-Dickson pairs in 
(A5). The conjugate quaternion corresponds to the hermitian conjugate q* = qR The quaternion norm squared is 
ijgjP = qq' = det q and the inverse quaternion q -1 = q' / det q. 
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We now superimpose the quaternionic structure on TV x TV matrices. Given four hermitian TV x TV matrices A, B,C, D 
we construct two non-lrermitian matrices X = A + iB and Y = C + iD and a 2TV x 2TV quaternionic matrix 

A 7~A ' T~\ ( A -)- iB C T %D \ ( X \ \ . , 

Q = A <g> CTo + B <g> 7 . 0-3 + c <g> ia 2 + D <g> 101 = ( _ c + iD J = ( _yt j^t ) - ( A >*0 • ( Al °) 

The upper row of blocks specify the whole matrix, so we can use the notation Q = ( X , Y ) is as in the Cayley-Dickson 
construction. As before, we define the functions F and S to extract the first and the second block of the pair 

X = F(Q) , Y = S(Q) . (All) 

Additionally we denote a trace in space of quaternionic matrices, which acts as block-trace operation in Pauli matrix 
representation, by TrbQ = (TrX, TrY). It projects quaternionic matrices to quaternions 

Tr bS = ( S) • (A12) 

Equipped with this matrix extension of quaternions one can naturally define pertinent objects to handle non-hermitian 
matrices including the quaternionic resolvent and the R transform. 


Appendix B 


We show here that Eq. ([5]) can be treated as a representation of the two-dimensional delta function. We have 


19 1 1 d p z — wj 19 2 1 \w\ 2 

IT dz Z + wj 7T dz \z\ 2 + |'U)| 2 IT dz \z\ 2 + \w\ 2 7T /, | 2 1 |2 


(Bl) 


The resulting function is a circularly symmetric bell-sliaped function located at the origin of the z-complex plane. 
The width of the peak at half maximum is of order |ui| and the height of order l/|w| 2 . For any non-zero value of w 
the integral of this function over the whole ^-complex plane is equal one 


/ 



1 . 


(B2) 


and it is independent of w. In the limit w —> 0 the volume under the peak stays constant while the width of the peak 
tends to zero. The whole function and its integral gets concentrated in a single point so it is the delta function. 

We note that also the second part of the inverse quaternion can be used to define the two-dimensional Dirac delta 

19^1 19 w 

7 T dw z + wj TT dw \z\ 2 + \w\ 2 

We see that when we fix the argument <f> of w = e l ^\w\ and take the limit |u>| —>• 0 + we obtain the delta function 
multiplied by a phase factor e 2l< ^ 2 ) (z). For <f> = 0 this expression reduces to the delta function. This representation is 
closer in spirit to the standard representation of the one-dimensional delta Q since the second part S for quaternions 
is like the imaginary part 3m for complex numbers. This representation is however more difficult to use in practice 
since one has to extract the dependence on w and calculate the derivative before one takes the limit. 


w 1 
w 7 r 


W 


■u; 


(B3) 


Appendix C 


One can show m that the expression in the brackets in Eq. (ITT]) behaves in the limit w —> 0 as a sum of delta 
functions located at eigenvalues of X 


lira —Tr H-Aw^HT 1 
N—too TV R 1 1 L 


1 

TV 


N 


^<5 (2) (z-A i) . 


(Cl) 
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We repeat this reasoning here. The two matrices Hl and Hr § are hermitian, therefore they have real eigenvalues. 
We assume that they are non-degenerate. Denote the corresponding eigenvectors by | L a ) and | R a ) respectively 


H L \L a ) = K a \L a ) 
Hr\Ro) = A q | R a ) 


(C2) 


The two matrices have the same eigenvalues. Indeed, multiplying the first equation on both sides by ( zl — X') and 
using the identity (zl — X^)Hl = Hr(z1 — X^) we get 


H R (zl - Xt)| L a ) = A Q (zl - A't)| L a ) , 


(C3) 


which means that \R a ) = 77(51 — X')\L a ) is an eigenvector of Hr to the same eigenvalue A Q . The coefficient 77 is to 
fix the norm of the vector (R a \R a ) = 1. Using the spectral decomposition of Hr and Hr one can cast the expression 
on the left hand side of Eq. (Cl) into the form 


1 

AT 


'117/ 


-l 


w 


1 ^V_ | oi j2|/r> ir \ |2 


= V £ £ 


N 


a=l 3=1 


M 2 i(iwr 

A„A« 


(C4) 


The eigenvalues A Q and eigenvectors \R a ) and \L a ) depend on X but also on w and z. In our notation this dependence 
is implicit. We will display it if needed. When z is not equal to an eigenvalue of X the eigenvalues A a ( z, w = 0) of 
H l and Hr © are strictly positive. As a consequence the expression on the right hand side of Eq. ( |C4[ ) behaves as 
|w| 2 for small w and vanishes for w = 0. The situation changes when z is equal to an eigenvalue of X because in this 
case the smallest eigenvalue A 0 (z = A j,w) = |w| 2 of Hl and Hr tends to zero for w —?• 0 and then the expression 


(C4) diverges as l/|u;| 


1 

N' 


Ir\ 


i 2 tj-1 „ i<i2oii 0 )r 
L m 


(C5) 


The divergence comes from the term a = [3 = 0 of the sum (C4). To summarize, in the limit w -3 0 the expression 
on the left hand side of Eq. © vanishes for all values of z except those equal to eigenvalues of A'. A more careful 
analysis of the behavior of the expression (C4) in the vicinity of eigenvalues of A shows that it behaves as delta 


functions located at those eigenvalues. In order to see this consider z close to an eigenvalue A j of A: such that the 
distance z — \ 3 is much smaller than the eigenvalue separation. One can apply the perturbation theory to determine 
the lowest eigenvalue of Hl and Hr Q. Up to the second order one finds 


A 0 « M 2 + \{R 0 \L 0 )\ 2 |z Xj 


(C6) 


Inserting A 0 to the diverging term (a = /3 = 0) of the sum (C4) and neglecting remaining terms, which are of order 
|w|° or |w| 2 , one finds for w —> 0 


^TrH^wfHZ 1 


\M 2 \(Ro\L 0 )\ 2 


N 


{\w\^ + \(Ro\Lo)\ 2 \z-X^ 


—5 (2) (z — Xj) . 


(C7) 


The same holds for z in the vicinity of any eigenvalue Aof X so one eventually arrives at Eq. (Cl). 


Appendix D 


In this appendix we recall diagrammatic derivation of Eqs. (151 and (28). We begin with hermitian matrices and 


towards the end we describe how to generalize the derivation to non-hermitian ones. 

Let us recall main ideas |441145| and introduce graphical notation. We consider hermitian random matrices defined 
by the probability measure 


1 


dn{H) oc DH exp —iVTrU (H) = DHexp-NTi ( -H 2 + —77 a + — H A + 


54 


(Dl) 


and apply Gaussian perturbation theory to calculate statistical averages (...) with respect to this measure. One does 
it by spliting the measure into the Gaussian part and the residual part V{H) = \H 2 + Vr(H), where Vr{H) = 
..., and treating the latter as a perturbation 


f H 3 + f H 4 


d»(H) = d^H) ( 1 + NTiVr(H) + i ( NTtV r (H )) 2 


(D2) 
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FIG. 7. Pictorial representation of the (a) propagator (D41 and (b,c) the first two vertices generated by the perturbative 
expansion ( |D2| l. In this appendix we use convention that external vertices are represented by open circles. In contrast, internal 
vertices (see below) are represented by filled circles. Internal vertices correspond to indices which are summed over. 


where dfj,o(H) oc exp — NTiH 2 /2 with a normalization ensuring that f dgo(H) = 1. Now the problem reduces to 
calculating averages of powers of H with respect to the Gaussian measure (.. .)o = / dgo(H) .... As follows from the 
Wick’s theorem, averages with respect to the Gaussian measure of higher order products of H can be replaced by 
products of the averages of second order products. For example for the fourth order the theorem tells us that 

{HabHcdHefHgh )o = (H ab H cd ) 0 (H e fH gh ) 0 + (H ab H e f) 0 {H cd H gh ) 0 + (H ab H gh ) 0 (H e fH cd } 0 . (D3) 

This observation underlies the idea of Feynman diagrams which are just graphical representation of products of 
factors ( H ab H cd } 0 . These factors are called propagators and are graphically represented as double lines between pairs 
of indices ab and cd as shown in Fig. [TJi. It is easy to find, by doing the Guassian integral, that the numerical value 
of the propagator for dgo(H) oc exp — NTtH 2 /2 is 


(HabH cd ) o — jyd ad S bc . 


(D4) 


This means that the only non-zero contributions of propagators are l/N for a equal to d and c to b. Feynman 
diagrams are obtained by connecting vertices (see Fig. p]),c) generated by perturbative expansion of the residual part 
Vr (if) (|D2[) with lines representing propagators (D4).^.n example of a diagram is shown in Fig. [S] A contribution 


of the diagram to the perturbative expansion is given by a product of contributions from propagators and vertices. 
Propagators contribute (D4) and vertices: Ng$/3 and Ng^/4 etc. The factor N comes from the prefactor before 
trace (D2). Let us concentrate on counting powers of N for a given diagram. Each propagator contributes a factor 
l/N each vertex contributes a factor N and each closed line, like the thick one drawn in Fig. |8j contributes a 

factor N. This is because the summation over internal indices of delta functions ( |D4| ) gives for a closed line a factor 
J2a = AT. Such a closed line can be viewed as a face of polyhedron made of edges (double lines) of the given 
diagram. The total power is thus N v + p - E where V is the number of vertices, F is the number of faces of the diagram 
viewed as polyhedron and E is the number of edges (represented us double lines). The combination V + F — E is the 
Euler characteristic of the polyhedron or equivalently of a two-dimensional surface on which the diagram can be drawn 
without edge-crossing. In the example shown in Fig. [8] we have V = 4, F = 4 and E = 6 and V + F — E = 2. This is 
equal to the Euler characteristic of sphere. Generally, V + F — E = 2 — 2h, where h is the genus of two-dimensional 
surface on which diagram can be drawn without crossings, h = 0 for sphere, h = 1 for torus, h = 2 for double torus, 
etc. Thus one can see that the leading contribution comes from planar diagrams, which can be drawn on sphere (or 
equivalently plane) without edge-crossing. Those which can be drawn on torus without edge-crossing are suppressed 
with l/N 2 factor, on double torus with l/N 4 etc. The diagrams which can be drawn on torus without edge-crossing 
can also be drawn on plane but then they have at least one edge-crossing. In Fig. [9]we show examples of lowest order 
diagrams which are planar ( h = 0) and non-planar (h = 1). One can generalize this classification to diagrams with 
external lines by considering two-dimensional surfaces with boundaries. 

The main conclusion is that in the limit N —> oo one can neglect contributions from non-planar diagrams which 
are suppressed with 1/N 2h powers as compared to the leading contribution coming from planar diagrams. In effect, 
in this limit one can restrict enumeration only to planar diagrams. This observation was first made in |43| and the 
technique of planar diagram enumeration was developed in |441 45! . 

So let us recall the idea of graphical enumeration of planar diagrams. Denote the n-point correlation function 


^i'23-2 ■ • ■ ^i n Jn) (D5) 

by a blob with n-external double legs as in Fig. |10f i. In perturbation theory the n-point correlation function is 
calculated as a sum over all planar Feynman diagrams with n external legs. An example of a diagram generated by 
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FIG. 8. Example of a planar Feynman diagram obtained from V = 4 cubic vertices and E = 6 propagators. It has F = 4 
closed lines, like e.g. the thick one surrounding the region 1. With each such closed line one can associate a face of polyhedron 
constructed from the edges of the diagram. If the diagram is drawn on a sphere, the region 4 can be viewed as a compact 
region or the fourth face of tetrahedron. 




the 4-point correlation function is shown in Fig. [UK One also defines an n-point connected correlation function 
which generates only connected diagrams with n external legs. We denote it by 






■ H, 


>» 


(D6) 


and in the graphical representation by a blob surrounded by a double contour (Fig. 10 a). An example of a diagram 
generated by the 4-point connected correlation function is shown in Fig. & Below we show how to use these 
functions in calculations of the Green’s function ([2]) . It is useful to consider a generalized resolvent which is an N x N 
matrix of the form 


G = 



(D7) 
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FIG. 10. Pictorial representation of (a) n-point correlation function fh n and (b) n-point connected correlation function K n . 





FIG. 11. Examples of diagrams contributing to 4-point correlation function. The diagram (b) contributes to 4-point connected 
correlation function. 


where the matrix Z is an arbitrary invertible constant hermitian matrix. Clearly, by taking the normalized trace 


G= lim —Tr G 

N—toe N 


(D8) 


and setting Z = zt one obtains the resolvent G(z) (|2j) . We will set Z = zt only at the end of the calculations, while 
they will be done for an arbitrary Z in order to keep track of the underlying algebraic structure. Writing the resolvent 
(D7| as a geometrical series one has 


G = ((Z- H)- 1 ) = Z~ l + ( Z~ X HZ~ X 


Z~ 1 HZ~ 1 HZ 


-l 


+ ... . 


(D9) 


The factors depending on Z can be written outside the brackets since Z is a constant matrix. What remains in the 
brackets are n-point correlation functions. For example, the contribution of the third term in the last equation to the 
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FIG. 12. Diagrammatic expansion of the resolvent, Eq. (|D9|) . Sum over internal indices is implicit. 



FIG. 13. Example of a diagram generated by resolvent. Its contribution is proportional to gz, gq and z 


matrix element G a f reads 

1^[Z- X HZ- X HZ,- X 


a/ 


= E (Za b 1 H bc Z- x H de Z~ x )= £ Z- b x Z- x Z: x (H b cH de ) . 


b.c.d.e 


(DIO) 


b.c.d.e 


Graphical representation of equation (D9) is shown in Fig. 12 The dashed lines correspond to elements of Z 1 and 


indic es are associated with the vertices at the endpoints of the dashed line. The internal vertices, like b 1 c, d and e in 
Eq. (DIO I are summed over, while the external ones a, / correspond to indices of the matrix element G a f ■ An example 
of a Feynman diagram generated by the resolvent G is shown in Fig. 13 This diagram is built out of two diagrams 
shown in Fig. [T4| sandwiched by dashed lines. Diagrams shown in Fig. [14] belong to a class of one-line-irreducible 
diagrams which are defined by the condition that they cannot be disconnected by removing a single dashed line. We 
denote a generating function of such diagrams by E which is an N x N matrix. Analogously, the indices of the E 
matrix correspond to the two external vertices and we define E (z) as a normalized trace of E: 


E = lim —Tr E. 

A^—foo N 


(Dll) 


The reason why it is convenient to introduce the generating function E is, that one can write down a set of diagram¬ 
matic equations relating E and G and derive from them a closed-form expression for the resolvent. These equations 
are often named after Dyson and Schwinger who invented a general class of such equations in field theory to sum 
contributions of different types of Feynman diagrams. 

In our case we will write down two Dyson-Schwinger equations to sum up all planar diagrams contributing to the 
resolvent G. The first Dyson-Schwinger equation exploits the fact that any diagram in G can be constructed as an 



FIG. 14. One-line-irreducible diagrams that are components of the diagram presented in Fig. [13] 
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-XXX . XXX .o + o. XXX . XXX . -XXX 


FIG. 15. Pictorial representation of the first Dyson-Schwinger equation (D12 I. Each diagram may be constructed by aligning 
some one-line-irreducible diagrams and connecting them alternately with dashed lines. 




FIG. 16. Pictorial representation of the second Dyson-Schwinger equation (D13). Each one-line-irreducible diagram may be 
constructed by aligning some diagrams next to each other and straddling them with a connected diagram. 


alternating chain of one-line-irreducible diagrams and dashed lines (compare Figs. 13 and 141. This leads to the 
graphical relation shown in Fig. 15 which is equivalent to the equation 


G = Z 


-l 


Z-^Z- 1 + Z~ L Y,Z-^Z 


7-1 


-1\ 


7-l\ 


7~1 


= (r-g) 


(D12) 


This is the first Dyson-Schwinger equation. There is another independent relation between G and S which exploits 
the fact that any one-line-irreducible diagram can be constructed by straddling G-diagrams by legs of a connected 
diagram as shown in Fig. 16 This graphical equation is equivalent to 


Xi a b — ( Kl)ab + £( ^2 )ac, dbGcd + ^ ' (^3 )ac,ef,dbG ce Gfd + ■ ■ ■ • 

cefd 


(D13) 


cd 


This is the second Dyson-Schwinger equation. The two Dyson-Schwinger equations simplify for Z = zt ( Z a b = zS a b) 
since in this case G = G(z)l and E = T,(z)l, where G(z) and E(z) are complex scalar functions. In effect, equations 
(D12| and (D13l reduce after taking the normalized trace to scalar equations 


G(z) = (z - E(*)) 


-1 


(D14) 


and 


E(^) — K\ + K 2 G(z) + ft3G 2 (z) + ... (D15) 

where 

Kn = lim ^«1WF*>> • (D16) 

at->oo Jy 

These coefficients are connected planar (or free) cumulants as they are defined as sums over all connected planar 
diagrams with n external legs. The generating function for planar cumulants is called R transform in free probability 
and it is defined as 

OO 

R(z ) = ki + k 2 z + K 3 Z 2 + ■ ■ ■ = y; K n z" _1 

n= 1 


(D17) 
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Using the R transform one can concisely write the second Dyson-Schwinger equation (D15) 


£(*) = R(G{z)). 

Inserting this into the first Dyson-Schwinger equation ( D14| ) one eventually obtains the formula 

G(2) = z-R(G{z)Y 


(D18) 


(D19) 


which was earlier mentioned in the main text (15). 


Generalization of these equations to non-hermitian matrices is straightforward. The Dyson-Schwinger equations 
(D12) and (D13) have the same form as for hermitian matrices but a slighly more complex index structure, since 


the N x N hermitian matrices are considered as quaternionic matrices (or equivalently 2N x 2TV matrices with a 
particular block structure) © and ( |18[ ). The additional index structure takes care of parts of quaternion (positions 
of blocks in the 2N x 2 N matrices). After taking the block-trace (A12)one gets a matrix equation for 2x2 matrices 
(quaternions) ([28]) instead of scalar equation (D19). 
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